#include <stdlib.h>
#include "dynamics.h"

int pack (int N, int status[N])
{
	int j, p;

	p = 0;
	for (j = 0; j < N; j++)
	{
		p = 2 * p + status[j];
	}

	return p;
}

void unpack (int p, int N, int status[N])
{
	int j;

	for (j = 0; j < N; j++)
	{
		status[N - 1 - j] = p % 2;
		p /= 2;
	}
}

void init_state (int N, int n, struct state state[n])
{
	int i;

	for (i = 0; i < n; i++)
	{
		state[i].status = calloc (N, sizeof (int));
		unpack (i, N, state[i].status);
	}
}

void reset_state (int n, struct state state[n])
{
	int i;

	for (i = 0; i < n; i++)
	{
		state[i].status0_first = -1;
		state[i].status0_last = -1;
		state[i].status1_sibling = -1;
		state[i].status1 = i;
		state[i].final = -1;
		state[i].flux = -1;
		state[i].flux2 = 0;
		state[i].dist1 = -1;
		state[i].dist2 = -1;
		state[i].basin = 0;
		state[i].overlap = 0.0;
		state[i].flag = -1;
	}                	
}                        	
                         	
/* calculate the next step for give state under the network using strong inhibition model */
int next_step (int N, int state[N], int network[N][N])
{
	int j, k, p, flag;
	int state2[N];

	for (j = 0; j < N; j++)
	{
		flag = 0;
		for (k = 0; k < N; k++)
		{
			if (k == j)
				continue;
			if (state[k] == 1 && network[j][k] == -1)
			{
				state2[j] = 0;
				flag = 1;
				break;
			}
		}

		if (flag == 1)
			continue;

		for (k = 0; k < N; k++)
		{
			if (k == j)
				continue;
			if (state[k] == 1 && network[j][k] == 1)
			{
				state2[j] = 1;
				flag = 1;
				break;
			}
		}

		if (flag == 1)
			continue;

		if (network[j][j] == -1)
			state2[j] = 0;
		else
			state2[j] = state[j];
	}

	p = state2[0];
	for (j = 1; j < N; j++)
	{
		p = 2 * p + state2[j];
	}

	return p;
}

/* get final state and distance to final state for a given state p[i] */
int find_final (int N, int i, void *p, int m)
{
	struct state *state;

	state = (struct state *) p;
	
	int n, r;

	n = state[i].status2;

	if (state[n].flag == -1)
	{
		state[i].flag = -2;
		r = find_final (N, n, state, m + 1);
		if (state[n].flag == -3) /* special treatment for the beginning of the cycle */
		{
			state[n].flag = 3;
			state[n].dist1 = 0;
			state[n].final = n;
			state[i].flag = 2;
			state[i].dist1 = 1;
			state[i].final = n;
			r = n;
		}
		else if (state[i].flag != -3 || m == 0)
		{
			state[i].flag = state[n].flag;
			if (state[n].flag == 3)
			{
				state[i].dist1 = 0;
				state[i].final = i;
			}
			else
			{
				state[i].dist1 = state[n].dist1 + 1;
				state[i].final = state[n].final;
			}
		}
	}
	else if (state[n].flag == -2) /* find a cycle, n is the start and i is the end  */
	{
		state[n].flag = -3;
		state[i].flag = 3;
		state[i].dist1 = 0;
		state[i].final = i;
		r = -1;
	}
	/* next state is a normal state or normal steady state */
	else if (state[n].flag == 0 || state[n].flag == 1)
	{
		state[i].flag = 0;
		state[i].dist1 = state[n].dist1 + 1;
		r = state[i].final = state[n].final;
	}
	/* next state is a cycle state or a state whose final state is a cycle state */
	else if (state[n].flag == 2 || state[n].flag == 3)
	{
		state[i].flag = 2;
		state[i].dist1 = state[n].dist1 + 1;
		r = state[i].final = state[n].final;
	}

	return r;
}

void generate_dynamics (int N, int network[N][N], int n, struct state state[n])
{
	int i, p;

	for (i = 0; i < n; i++)
	{
		p = next_step(N, state[i].status, network);
		state[i].status2 = p;
        
		if (state[i].status1 == state[i].status2)
		{
			state[i].flag = 1;
			state[i].final = i;
			state[i].dist1 = 0;
		}
		else
		{
			if (state[p].status0_first == -1)
			{
				state[p].status0_first = i;
			}
			else
			{
				state[state[p].status0_last].status1_sibling = i;
			}
			state[p].status0_last = i;
		}
	}

	for (i = 0; i < n; i++)
	{
		if (state[i].final == -1)
		{
			find_final (N, i, state, 0);
		}
		state[state[i].final].basin ++;
	}
}

/* calculate the total flux for a given state p[i] */
int calc_flux2 (int N, int i, void *p)
{
	struct state *state;

	state = (struct state *) p;

	int r, flux;

	r = state[i].status0_first;
	flux = 0;
	
	while (r != -1)
	{
		if (state[r].flux == -1)
			flux += calc_flux2 (N, r, state) + 1;
		else
			flux += state[r].flux + 1;

		r = state[r].status1_sibling;
	}

	return flux;
}

void calc_flux (int N, int n, struct state state[n])
{
	int i;

	for (i = 0; i < n; i++)
	{
		if (state[i].flag == 0 || state[i].flag == 2)
		{
			state[i].flux = calc_flux2 (N, i, state);
			state[i].flux1 = state[i].flux + 1;
		}
		else if (state[i].flag == 1 || state[i].flag == 3)
		{
			state[i].flux = state[i].basin - 1;
			state[i].flux1 = state[i].basin;
		}
	}
}
